library(ggplot2)
library(countrycode)
library(rworldmap)
library(tidyr)
library(tidyverse)
library(cutpointr)
library(devtools)
library(GGally)
library(dplyr)
library(systemfit)

#dummies package is not supported with newer versions of R. This downloads the version used for the paper 
install_version("dummies", version ="1.5.6", repos = "http://cran.us.r-project.org")
library(dummies)

options(stringsAsFactors = FALSE)
options(scipen=999)

#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  


#load final data set and worker's rights source data 
final_dat <- read.csv("~/Desktop/final_dat.csv")
dat<- final_dat 

workerr.latent <- read.csv("~/Desktop/workerr_latent.csv")

#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  

#### Figure 1 ###### 
#Subset to Countries of Interest
#U.S., U.K., Switzerland, Sweden, Spain, Portugal, Norway, Netherlands, Japan, Italy, Ireland, Greece, Germany, France, Finland, Denmark, Canada, Belgium, Austria, Australia 
#countrycode package makes it easy to merge and switch between ccode, name, and iso3 

#create a temp dataframe with subset of countries for Figure 1 and Figure 2, add country names for plots 
temp <- dat[dat$ccode %in% c(2, 740, 255, 220, 325, 200, 20, 305, 211, 390, 375, 350, 205, 210, 385, 235, 230, 380, 225, 900), ] 
temp$country <- countrycode(temp$ccode, "cown", "country.name")

#Boxplot of PIR
ggplot(temp,aes(x=country, y=latentmean))+
  geom_boxplot() + 
  theme_bw() + 
  labs(x= "Country", y="Latent PIR Score") + 
  coord_flip()

#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  

#### Figure 2 ###### 
#Subset to 2002 
temp <- temp[temp$year == 2002, ]

#create dataframe for PIR that only contains the necessary columns - country name and PIR value, adds a description column for PIR value
foo<- data.frame(temp$country, temp$latentmean, variable="PIR")
colnames(foo) <- c("country", "value", "variable") 
foo$country <- factor(foo$country, levels = foo$country[order(foo$value)])

#create dataframe for worker rights that only contains necessary columns and adds a description column for workers rights 
bar<- data.frame(temp$country, temp$workerr.mean, variable="Workers' Rights") 
colnames(bar) <- c("country", "value", "variable") 

#combine foo and bar to make one long dataframe 
plot_dat <- rbind(foo, bar) 

#plot side by side comparison 
ggplot(plot_dat, aes(x = value, y = country)) +
  geom_point(alpha=1) + 
  facet_wrap(~ variable, scales = "free_x") + 
  labs(x = "", y= "Country") + 
  theme_bw()
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  


#### Figure 3 ###### 
#create map 

#create a temporary working dataframe - dat_map
dat_map <- final_dat

#subset to 1995 
dat_map <- dat_map[dat_map$year==1995,]
dat_map <- dat_map[dat_map$ccode!=315, ] # no data for 315 (Czechoslovakia), creates an error later if kept in

#rworldmap likes iso3, so add that variable from cowcodes 
dat_map$iso3 <- countrycode(dat_map$ccode, origin="cown", destination = "iso3c")

#keep just access points and iso3 data 
dat_map <- dat_map[,c("access_points","iso3")]

#remove missing data 
dat_map<- dat_map[complete.cases(dat_map[,c("access_points", "iso3")]), ]

#creating data for rworldmap package to read 
joinData<- joinCountryData2Map(dat_map,
                               joinCode="ISO3", 
                               nameJoinColumn = "iso3")

#create plot 
theMap <- mapCountryData(joinData, nameColumnToPlot = "access_points", addLegend = TRUE, catMethod = "fixedWidth", colourPalette = "white2Black", mapTitle="Access Points by Country")

#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  

###########Table 1###########
#Correlation between component parts and the workers' latent variable 

worker <- workerr_latent[,-c(1,2,3,36,38)]
cor(worker$minimum.wage, worker$workerr.mean, use = "complete.obs")
cor(worker$association.govt, worker$workerr.mean, use = "complete.obs")
cor(worker$association.mrkt,  worker$workerr.mean, use = "complete.obs")
cor(worker$colbargain.govt,  worker$workerr.mean, use = "complete.obs")
cor(worker$colbargain.mrkt, worker$workerr.mean, use = "complete.obs")
cor(worker$strike.govt,  worker$workerr.mean, use = "complete.obs")
cor(worker$strike.mrkt,  worker$workerr.mean, use = "complete.obs")
cor(worker$ciri_worker,  worker$workerr.mean, use = "complete.obs")
cor(worker$LaborRightsPos,  worker$workerr.mean, use = "complete.obs")
cor(worker$LawPos,  worker$workerr.mean, use = "complete.obs")
cor(worker$PracticePos,  worker$workerr.mean, use = "complete.obs")
cor(worker$uri_l,  worker$workerr.mean, use = "complete.obs")
cor(worker$uri_p,  worker$workerr.mean, use = "complete.obs")
cor(worker$sri_l,  worker$workerr.mean, use = "complete.obs")
cor(worker$sri_p, worker$workerr.mean, use = "complete.obs")


#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  
#----------------------------------------------------------------------------------------------------------------------------------------------------------------#  


########## Table 2################
dat<- final_dat 

#create dummy variables for fixed effects and merge with data
country.dummies <- as.data.frame(dummy(dat$ccode)) #produces a warning due using T/F value for contrasts arguments instead of a matrix of contrasts. Thi of R no longer ignores and throws a warning https://stackoverflow.com/questions/56637183/warning-message-dummy-from-dummies-package
colnames(country.dummies) <- c("ccode2", "ccode20" ,"ccode31" ,"ccode40" ,"ccode41" ,"ccode42" ,"ccode51" ,"ccode52" ,"ccode53" ,"ccode54" ,"ccode55" ,"ccode56" ,"ccode57" ,"ccode58" ,"ccode60" ,"ccode70" ,"ccode80" ,"ccode90" ,"ccode91" ,"ccode92" ,"ccode93" ,"ccode94" ,"ccode95" ,"ccode100" ,"ccode101" ,"ccode110" ,"ccode115" ,"ccode130" ,"ccode135" ,"ccode140" ,"ccode145" ,"ccode155" ,"ccode160" ,"ccode165" ,"ccode200","ccode205" ,"ccode210" ,"ccode211" ,"ccode212" ,"ccode220" ,"ccode223" ,"ccode225" ,"ccode230" ,"ccode232" ,"ccode235" ,"ccode255" ,"ccode290" ,"ccode305" ,"ccode310" ,"ccode315" ,"ccode316" ,"ccode317" ,"ccode325" ,"ccode331" ,"ccode338" ,"ccode339" ,"ccode343" ,"ccode344" ,"ccode349" ,"ccode350" ,"ccode352" ,"ccode355" ,"ccode359" ,"ccode360" ,"ccode365" ,"ccode366" ,"ccode367" ,"ccode368" ,"ccode369" ,"ccode371" ,"ccode375" ,"ccode380" ,"ccode385" ,"ccode390" ,"ccode395" ,"ccode402" ,"ccode403" ,"ccode404" ,"ccode432" ,"ccode433" ,"ccode434" ,"ccode436" ,"ccode437" ,"ccode451" ,"ccode452" ,"ccode475" ,"ccode482" ,"ccode484" ,"ccode500" ,"ccode501" ,"ccode516" ,"ccode520" ,"ccode551" ,"ccode553" ,"ccode560" ,"ccode565" ,"ccode570" ,"ccode580" ,"ccode581" ,"ccode590" ,"ccode625" ,"ccode640" ,"ccode660" ,"ccode666" ,"ccode712" ,"ccode713" ,"ccode732" ,"ccode740" ,"ccode750" ,"ccode770" ,"ccode771" ,"ccode775" ,"ccode780" ,"ccode790" ,"ccode800" ,"ccode812" ,"ccode840" ,"ccode850" ,"ccode900" ,"ccode910" ,"ccode920" ,"ccode935" ,"ccode940" ,"ccode946" ,"ccode970" ,"ccode983" ,"ccode986" ,"ccode987")  
dat <- cbind(dat, country.dummies) 

#set model equations 
latentmean.model<- latentmean ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity #model 1 
workerr.model <-  workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil #model 1
latentmean.model.fe <-latentmean ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.log.oil + l.conflict + l.cat_ratify + l.polity  +   ccode41  +  ccode42  +  ccode51  + ccode90 + ccode91 + ccode92 + ccode93 + ccode94 + ccode95 + ccode100 + ccode101 +   ccode130 + ccode140 + ccode145  + ccode160  + ccode339 +  ccode355 + ccode360  +  ccode432  + ccode434 + ccode436 +  ccode451  + ccode482 + ccode551 + ccode553 + ccode560  + ccode590 +  ccode640 +  ccode712 +  ccode750 + ccode771 + ccode780 + ccode790 + ccode800 + ccode840 #model 2 (with fixed effects, removing the ones that cause separation issues)
workerr.model.fe <- workerr.mean2  ~ l.access_points + l.log.realgdppccurrent + l.log.pop + l.ythblgap + l.trade_WDI + l.fdiflows_UNCTAD  + l.conflict  +  l.under_BCG + l.cescr_ratify + l.polity + l.log.oil   +  ccode41  +  ccode42  +  ccode51  + ccode90 + ccode91 + ccode92 + ccode93 + ccode94 + ccode95 + ccode100 + ccode101 +  ccode130 + ccode140 + ccode145  + ccode160  + ccode339 +  ccode355 + ccode360  +  ccode432  + ccode434 + ccode436 +  ccode451   + ccode482 + ccode551 + ccode553 + ccode560+ ccode590 +  ccode640 +  ccode712 +  ccode750 + ccode771 + ccode780 + ccode790 + ccode800 + ccode840 #model 2 (with fixed effects, removing the ones that cause separation issues)

#create complete cases of observations 
dat1 <- dat[complete.cases(dat[,c("latentmean", "l.access_points", "l.log.realgdppccurrent", "l.log.pop", "l.ythblgap", "l.trade_WDI", "l.log.oil","l.conflict", "l.cat_ratify", "l.polity", "l.fdiflows_UNCTAD", "l.lji_LS", "l.cescr_ratify", "l.under_BCG", "workerr.mean")]), ]

#fit SUR (model 1)
fitsur1 <- systemfit(list(surfit.latentmean1=latentmean.model, surfit.workerr1=workerr.model), method="SUR", data=dat1)
summary(fitsur1)

#fit SUR (model 2)
fitsur2 <- systemfit(list(surfit.latentmean1=latentmean.model.fe , surfit.workerr1=workerr.model.fe), method="SUR", data=dat1)
summary(fitsur2)


